Dynamical system analyser

ABSTRACT

A dynamical system analyser (10) incorporates a computer (22) to perform a singular value decomposition of a time series of signals from a nonlinear (possibly chaotic) dynamical system (14). Relatively low-noise singular vectors from the decomposition are loaded into a finite impulse response filter (34). The time series is formed into Takens&#39; vectors each of which is projected onto each of the singular vectors by the filter (34). Each Takens&#39; vector thereby provides the co-ordinates of a respective point on a trajectory of the system (14) in a phase space. A heuristic processor (44) is used to transform delayed co-ordinates by QR decomposition and least squares fitting so that they are fitted to non-delayed co-ordinates. The heuristic processor (44) generates a mathematical model to implement this transformation, which predicts future system states on the basis of respective current states. A trial system is employed to generate like co-ordinates for transforation in the heuristic processor (44). This produces estimates of the trial system&#39;s future states predicted from the comparison system&#39;s model. Alternatively, divergences between such estimates and actual behavior may be obtained. As a further alternative, mathematical models derived by the analyser (10) from different dynamical systems may be compared.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates to a dynamical system analyser, and moreparticularly to such a device applicable to analysis of dynamicalsystems which might be nonlinear or chaotic in the mathematical sense.

2. Discussion of Prior Art

Investigation of nonlinear dynamical systems presents particulardifficulties because traditional electronic spectral analysis tools arenot appropriate. As an example, consider a chaotic nonlinear dynamicalsystem arranged to produce an electronic signal characteristic of itsbehaviour. It is known from dynamical systems research that such asignal can have a substantial continuous component to its spectrum.Filtering in either time or frequency to improve signal to noise ratioscan distort or alter the system dynamics as perceived from a filteredsignal. This is exemplified by the papers of Badii et al in Phys. Rev.Lett. 60 (1988) p.979 and Milschske et al. in Phys. Rev. A37 (1988) p.4518. Moreover, power spectrum analysis is insufficient to characterisethe dynamics of a system when data exhibits deterministic chaosresulting in a continuous spectrum (as opposed to a set of harmonicswith a simple relationship). This is discussed by F C Moon in a workentitled "Chaotic Vibrations", Wiley-Interscience. Spectrally selectivesignal processing of this kind, intended to simplify analysis, canrender the apparent system behaviour more complex.

The behaviour of a dynamical system is commonly represented by a curvein a multi-dimensional phase space. Successive points on the curveindicate the evolution of the system as a function of time. The line isreferred to as a "trajectory", and the region of phase space to which itis confined after an arbitrary time is called an "attractor". If thedynamical system is chaotic, the region is called a "strange attractor".The attractor is the phase space region in which the system behaviour islocalised, and to which it can be said to be "attracted".

There is a need for a device which is applicable to the analysis ofnonlinear systems, since traditional techniques merely complicate theapparent behaviour instead of reducing it to identifiable components.Traditional techniques Such as Fourier spectral analysis are applicableonly to systems which can be modelled as the linear superposition of alimited number of harmonic modes.

There is a particular need for a device which can detect a change fromnormal system behaviour to chaotic behaviour. This is of particularimportance in the field of aero-engines and other mechanical motivepower sources. Such power sources are characterised by highly regularoscillatory and/or rotational behaviour in normal operation, but developirregular characteristics before failing catastrophically. A devicecapable of detecting onset or wear-related irregular behaviour wouldprovide a means for anticipating and avoiding catastrophic failure byprior shut-down and repair. It might also provide for reduced servicingcosts, since maintenance could be deferred until system behaviourwarranted it.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide a dynamical systemanalyser.

The present invention provides a dynamical system analyser characterisedin that it includes:

(1) means for generating a sequence of sets of phase space coordinatesfrom a time series of signals from a dynamical system, each coordinateset being projections of a respective Takens' vector on to a set ofsingular vectors obtained in a singular value decomposition of the saidtime series of signals or of another such series, and

(2) heuristic processing means arranged to carry out a transformation ofthe said sequence to produce a fit to reference data and to create amathematical model related to that transformation.

The invention provides the advantage that it produces a mathematicalmodel of a dynamical system, and is applicable to the analysis ofnon-linear systems. As will be described later, a mathematical modelcreated by the heuristic processing means may be employed to generatepredictions of future dynamical system behaviour on the basis ofpreceding behaviour. Alternatively, mathematical models derived fromdifferent dynamical systems may be compared with one another for thepurposes of similarity or change analysis.

In a preferred embodiment, the invention provides a dynamical systemanalyser including:

(1) means for deriving a time series of signals from a dynamical system,

(2) means for generating from the time series a set of singular vectorseach associated with a respective singular value which is greater than anoise level of the dynamical system, the set defining a phase space andcorresponding to a subset of a set of vectors from a singular valuedecomposition of the time series,

(3) means for transforming the time series into a sequence of sets ofphase space co-ordinates, and

(4) heuristic processing means arranged to operate in a training mode tocarry out a transformation in which the said sequence undergoes QRdecomposition (as herein defined) and least squares fitting to referencedata in order to generate a mathematical model related to thattransformation.

The heuristic processing means may also be arranged to operate in a testmode to transform a further sequence of phase space co-ordinates, and toprovide predictive output values on the basis of the mathematical model.

In this embodiment, the invention provides the advantage of beingtrainable in accordance with reference data to generate a mathematicalmodel. This enables reference data to be selected and also to bereplaced if retraining is required.

The heuristic processing means may be arranged to employ in trainingmode reference data which is either one of or both of the time seriesand the sequence of sets of phase space co-ordinates derived therefrom.The reference data is input to heuristic processing means with arelative lead time interval over input of the said sequence for QRdecomposition and least square fitting; the lead time interval isarranged to be not greater than the inherent predictable timescale ofthe dynamical system as indicated by its largest Lyapunov component. Inthis case, the mathematical model which is generated in training mode isthat required to extrapolate the phase space trajectory and/or the timeseries of the dynamical system ahead in time to the extent of the leadtime interval. In test mode, the output values correspond to predictionsof future behaviour of a trial system which produced the furthersequence of phase space co-orinates previously referred to. These valuesmay be the phase space co-orinates and/or the time series of that trialsystem predicted on the basis of the mathematical model obtained intraining mode.

Alternatively, or in addition, the heuristic processing means may bearranged to compare the mathematical model's predictions of the trialsystem's phase space co-orinates and/or time series with actual valuesof these parameters measured later in time. In this case, error valuesarise between each prediction and a respective measured parameter. Alarge error value indicates significant departure from a predictedextrapolation of a phase space trajectory. In a mechanical system, thismight be associated with impending catastrophic failure or otherunwelcome behaviourial change.

The generating means may be a digital computer programmed to carry out asingular value decomposition of the dynamical system time series. Thismay be a once and for all computation, in which case its results may bestored in a memory for use when required.

The dynamical system may be a system (such as a mechanical or electronicdevice) in an initial functionally acceptable state. A trial system tobe compared with the dynamical system may be the same system afterdegradation has occurred in service over a period. This would enable themaintenance requirements of large machines (for example) to be assessedon the basis of behaviourial change.

The transforming means may be a finite impulse response (FIR) filteringdevice arranged to input the time series in the form of Takens' vectors(as hereinafter defined), and to compute the projections of each Takens'vector on to each of the set of singular vectors obtained by singularvalue decomposition to provide a respective set of co-ordinates in aphase space for each such Takens' vector.

The heuristic processing means is preferably substantially as describedin an International Patent Application No PCT/GB90/00142 published underthe Patent Co-operation Treaty as No WO 90/09643.

The dynamical system analyser of the invention may be employed with theheuristic processing means operating in training mode to providemathematical models of respective dynamical systems. The models may thenbe compared to provide an indication of the degree of similarity betweenthe systems. Singular vectors obtained from one dynamical system's timeseries may be employed in the production of phase space coordinates foreach of the systems' time series. Alternatively, each system's timeseries may be used in both singular value decomposition andtransformation by the heuristic processing means.

BRIEF DESCRIPTION OF THE DRAWINGS

In order that the invention might be more fully understood, anembodiment thereof will now be described, with reference to theaccompanying drawings, in which:

FIG. 1(a)-(b) is a schematic block diagram of a dynamical systemanalyser of the invention in the form of a prediction device;

FIG. 1(b) a shows a part of FIG. 1(a) in more detail;

FIG. 2 is a graph of a time series of signals for input to the

FIG. 1(a) device;

FIG. 3 is a graph of an autocorrelation function of the FIG. 2 timeseries;

FIGS. 4(a) and 4(b) provides the results of a singular valuedecomposition of the FIG. 2 time series and indicates a system noiselevel above digital quantisation noise;

FIG. 5(a) and 5(b) are is equivalent to FIG. 4 except that the relevantsystem noise level is below digital quantisation noise; FIG. 6schematically illustrates an alternative embodiment of the inventionarranged for determining system predictability; and

FIG. 7 is a graph showing results of the kind which might be obtainablefrom the FIG. 6 embodiment.

DETAILED DISCUSSION OF PREFERRED EMBODIMENTS

Referring to FIGS. 1(a) and 1(b), there is shown a schematic blockdiagram of a dynamical system analyser device of the invention deviceindicated generally by 10.

The device 10 incorporates a sensor transducer 12 mounted on anon-linear system 14 and deriving an analogue electrical signal A(t)characteristic of that system's behaviour. Output signals from thesensor 12 are fed to an analogue to digital converter (ADC) 16. Thelatter has an input S/H arranged to sample and hold the signal A(t) atsuccessive times t₁ . . . t_(i) . . . t_(M) with a constant spacing τ.Here t_(i) is given by:

    t.sub.i =t.sub.1 +(i-1)τ                               (1)

The analogue signal A(t) and sampling times t₁ etc are shown in a graph18. Discrete samples of the signal form a time series A(t_(i)), wherei=1 to M. These are digitised by the ADC 16 to produce a digital timeseries V(t_(i)) where V(t_(i)) corresponds to A(t_(i)) in each case fori=1 to M.

The ADC 16 operates under the control of a sampling clock 20, whichsupplies timing pulses at times τ apart. The magnitude of τ is set by acomputer 22, which supplies a digital word to the clock 20 indicatingthe required sampling interval.

The digital time series with general term V(t_(i)) is output from theADC 16 to the computer 22 in bit-parallel form via a computer input bus24. The computer 22 stores this time series on a hard disc, andcalculates the mean and variance of the time series terms V(t_(i)) toV(t_(M)). It produces a normalised version of the time series, and thenprocesses it to provide a singular value decomposition. Thedecomposition will be described in more detail later. It yields a set ofsingular vectors f₁ to f_(d), each of which is a mathematical functionexpressed as a series of n successive values. Some of the vectors areillustrated schematically in boxes such as 28, and their correspondingsingular value spectrum is indicated schematically at 29.

The computer 22 loads the singular vectors f₁ to f_(d) into a randomaccess memory (RAM) 30 via a first output bus 32. The RAM 30 hasrespective portions 30₁ to 30_(d) for the vectors f₁ to f_(d). Eachportion such as 301 stores n values per vector, eg f₁₁ to f_(1n) for f₁.The RAM 30 is incorporated in a multi-dimensional finite impulseresponse (FIR) filter 34, which also incorporates a chain 36 of nmulti-bit delay latches 36₁ to 36_(n). Parts of the RAM 30 and latchchain 36 which are not illustrated explicitly are indicated by dottedlines. The FIR filter 34 is of generally known and commerciallyavailable kind. One such incorporating a RAM 30 and latch chain 36 wasoffered for sale in May 1985 by Calmes Systems Inc., a Canadiancorporation. The FIR filter 34 may alternatively be based on devices asset out in U.S. Pat Nos. 4,885,715 and 4,833,635 referred to ascorrelators and convolvers.

The latch chain 36 receives input of the stored digital time seriesvalues such as V(t_(i)) from the computer 22 via a spur 38a of a secondoutput bus 38. The bus 38 also supplies a clock signal to the FIR filter34. This provides for successive-values V(t_(i)) etc in the digital timeseries to be clocked into a final latch 36_(n) of the chain 36. Eachvalue is a multibit parallel digital word, and its bits are inputsynchronously. After input, each time series value is clocked onwardsdown the chain towards a first latch 36₁ at the rate of one latch perclock cycle. After n clock cycles, the digital time series portion V(t₁)to V(t_(n)) is stored on latches 36₁ to 36_(n) respectively. This isknown in electronics as a delay register prefill operation. Afterprefill, the computer 22 supplies a clock signal actuation to the FIRfilter 34 to initiate its filtering operation. The operation of thefilter 34 is to multiply each time series value by a respective singularvector horizontally opposite in FIG. 1(a) to form products such asV(t_(i))f₁₁ through to V(t_(i+n-1))f_(1n). The products formed on eachclock cycle are added to provide a respective total for each singularvector.

Defining the first filter clock cycle as that during which the digitaltime series signal or element V(t₁) became loaded into the FIR filter34, then the filter outputs on the ith filter clock cycle co-ordinatesg₁ (i) to g_(d) (i) given by: ##EQU1## Equations (2.1) to (2.d) indicatethat the ith digital time series value V(t_(i)) is used by the FIRfilter 34 to form a set of co-ordinates g₁ (i) to gd(i). Theseco-ordinates are produced over the time series subset of n successivevalues V(t_(i)) to V(t_(i+n-1)), and represent the position at timet_(i) of a vector representing the nonlinear system in a phase spacehaving d dimensions defined by respective basis vectors f₁ to f_(d). Thegeneral co-ordinate gk(i) is the projection of the time series subsetV(t_(i)) to V(ti+n-₁) on to the general basis vector (singular vector)f_(k). The subset of n values V(t_(i)) to V(ti+_(n-1)) is referred to asa Takens' vector, after F Takens, "Detecting Strange Attractors inTurbulence", Lecture Notes in Mathematics, Ed. D. A. Rand and L. S.Young, Springer, Berlin 1981 p.366. The ith Takens' vector has an lthelement T₁ (i) defined by:

    T.sub.1 (i)=V(t.sub.i+1-1)                                 (3)

Varying the element number 1 from 1 to n generates V(t_(i)) toV(t_(i+n-1)). The general co-ordinate g_(k) (i) produced on the ithfilter clock cycle may be considered as the correlation of the ithTakens' vector with the sequence of values f_(k1) to f_(kn) of the kthsingular vector. Mathematically it is the projection of the ith Takens'vector on to the kth singular vector.

The co-ordinates g₁ (i) to g_(d) (i) are output in parallel and insynchronism from the FIR filter 34 on a filter bus 40. Successiveco-ordinate sets g₁ (i+1) to g_(d) (i+1), g₁ (i+2) to gd(i+2) etc areoutput on successive clock cycles. The filter bus 40 has a first spur40a connected to a multipole, double throw switch 42, and thence to ananswer input AI of an heuristic processor indicated generally by 44.Here the legend "HEURISTIC PROCESSOR" indicates elements AI, DI and 50to 68b described below. The throws of the switch 42 connect the answerinput AI either to the first bus spur 40a or to a zero signal (notshown). A second spur 40b of the filter bus 40 is connected to a delayunit 46, which delays each of the coefficients g₁ (i) to g_(d) (i) bylike time intervals of pτ. Here p is an integer and τ is the length of aclock cycle, so the delay is the same integral number of clock cyclesfor each of the co-ordinates g₁ (i) to g_(d) (i).

The delay unit 46 is connected by a delay output bus 48 to a data inputbus DI of the heuristic processor 44. A second spur 38b of the secondcomputer output bus 38 is connected via the switch 42 to a subsidiaryanswer input bus SAI of the heuristic processor 44. An input interfaceof the processor 44 is indicated by a chain line 50.

The heuristic processor 44 is of known kind. A detailed description ofits construction and mode of operation is given in InternationalApplication No PCT/GB90/00142 published under the Patent Co-operationTreaty (PCT) on 23 August 1990 with an International Publication No WO90/09643. Its features will be discussed briefly.

The heuristic processor 44 includes a digital arithmetic unit 52. Thelatter incorporates arithmetic circuits P₁₁ to P_(dr) shown asrectangles. These are arranged in rows P₁₁ to P_(1r),. . . P_(d1) toP_(dr) and columns P₁₁ to P_(d1), . . . P_(1r) to Pdr- The parameter dis the dimensionality of the co-ordinates g₁ (i) etc as before. Theparameter r is the number of dimensions of the space into which datainput to the unit 52 is transformed. The unit 52 also incorporatesmemories M₁ . . . M_(r) shown as diamonds. These are each disposed toreceive output signals from a respective column of arithmetic units,such as P₁₁ to P_(d1) for memory M₁ for example. Devices of the kind Por M which are not illustrated explicitly are denoted by chain lines.

Output from the digital arithmetic unit passes to a triangular systolicarray 54 which is also part of the heuristic processor 44. Thetriangular array 54 is shown vertically foreshortened for illustrationalconvenience. It incorporates boundary cells B₁₁ to Brr indicated bycircles, together with internal cells, I₁₂ to I_(r-1),r indicated bysquares. Structure not illustrated explicitly is denoted by chain lines.Each boundary cell (referred to generally as B) of the triangular array54 is arranged to calculate rotation parameters from data input fromabove, to update a respective stored matrix element and to output theparameters laterally to the right. Each internal cell (referred togenerally as I) of the triangular array 54 is arranged to receiverotation parameters from the left, to apply them to data received fromabove, to generate an output below, to update a respective stored matrixelement, and to pass on the rotation parameters to the right. Thesefunctions take place on each clock cycle, the heuristic processor 44receiving a clock signal from the computer 22 via a line 56. This clocksignal is relayed to each element of the heuristic processor 44.

The boundary and internal cells B and I are transputers programmed tocarry out arithmetic functions. They are designed to implement arotation algorithm. A matrix of transformed data is output from thedigital arithmetic unit 52 at the rate of one matrix row per clockcycle. This is input to the triangular array 54 at a like rate, andelements of the rows are progressively eliminated by a series ofrotations in accordance with the algorithm at the rate of one elementper row. The result is that the matrix of transformed data from the unit52 undergoes rotations which convert it to a triangular matrix. Thereare a number of rotation algorithms which the triangular matrix 54 arraymight employ. One such employs pairs of rotation parameters consistingof the sines and cosines of respective rotation angles. The boundary andinternal cells B and I store and update matrix elements of thetriangular matrix. Other computationally more convenient algorithmsproduce rotation parameters and stored elements mathematically relatedto such sines, cosines and matrix elements. In particular, in onealgorithm, the boundary cells B store a matrix of a first kind and theinternal cells I store a matrix of a second kind intended formultiplying the first kind matrix to generate the triangular matrix.

The process of triangularising a matrix by a series of plane rotationsis known as QR decomposition, Q being a matrix of rotation parametersand R the triangular matrix. Triangular array processing functions forimplementing rotation algorithms are well known in the prior art, suchas the international application referred to above, and will not bedescribed further.

The rotation parameters pass along rows of cells such as uppermost rowB₁₁ to I_(1r). They thereafter pass into a rectangular array 58 ofinternal cells I₁,r+1 to I_(r-1) and multiplier cells M_(r+1),r+1 toM_(r+1),r+d+1 indicated by hexagons. The boundary cells B₁₁ to Brr havethe additional function of cumulatively multiplying their cosinerotation parameters (or cosine-like parameters in other algorithms). Forexample, when the first boundary cell B₁₁ generates sine and cosineparameters on a clock cycle and outputs them to internal cell I₁₂ to theright, it also outputs the cosine parameter diagonally below and to theright to the second row, second column boundary cell (not shown, wouldbe designated B₂₂). This cell multiplies by its cosine parametergenerated two cycles later, there being a two clock cycle delay betweena diagonal output from one boundary cell and a like output from itsdiagonal neighbour incorporating the first such cell's Output. Thisdiagonal boundary cell function generates products arising fromcumulatively multiplied cosine parameters, each product generated from arespective input row of the matrix output from the digital arithmeticunit 52. Each product is fed to the row of multiplier cells M_(r+1) toM_(r+1),r+d+1, and passes along it at the rate of one multiplier cellper clock cycle. Similarly, rotation parameters passing along rows ofthe triangular array 54 are fed into respective rows of the rectangulararray 58. These parameters are employed by the rectangular array' sinternal cells I₁,r+1 etc to rotate data received from above, to updaterespective stored elements, provide output below and pass the parametersto the right; ie the rectangular array's internal cells operate inprecisely the same way as those of-the triangular array 54.

The multiplier cells (referred to generally as M) have the processingfunction of receiving data from above, multiplying it by an input fromthe left, outputting the product below and passing on to the right theleft hand input. The product is the scalar least squares residual errorobtained in fitting the transformed data from the arithmetic unit 52 todata passing down the column surmounting the relevant multiplier cell ineach case. Sets of residuals output by the rectangular array 58 formrespective least squares residual vectors.

The rectangular array 58 is notionally partitioned into a multi-columnportion 58a and a single column portion 58b.

The triangular array 54 and the rectangular array 58 receive differentinput data. The former receives data from the input DI which hasundergone transformation in the digital arithmetic unit 52, as has beensaid. The rectangular array receives input data from the answer input AIvia chains of latches (not illustrated explicitly). These are indicatedas an L-shaped delay device 60 incorporated in the heuristic processor44; the device 60 is notionally partitioned into a multichannel portion60a and a single channel portion 60b connected respectively to portions58a and 58b of the rectangular array 58.

The data input bus DI is connected to the arithmetic unit 52 via atemporal skewing device indicated by a triangle 62. This device is atriangular array of clocked latches (not shown). It introduces a timedelay which varies in magnitude across the input bus DI. The delay iszero at the upper edge of the input bus DI, and increases by one clockinterval τ per bus segment (not illustrated) to (d-1)τ at the lower edgeof the bus. Similarly, a second temporal skewing device 64 introduces alike skew delay variation at 64a in the answer input bus AI togetherwith a delay of dτ at 64b in the subsidiary answer input bus SAI. Outputfrom the rectangular array 58 passes via a temporal deskewing device 66similar to the device 62. The former introduces a rectangular arrayoutput delay varying from dτ for the column of cells I₁,r+1 toM_(r+1),r+1 to zero at the rightmost column of cells II,r+d+₁ toM_(r+1),r+d+1. The delay reduces by one clock cycle interval τ percolumn.

The purpose of the input delay devices 62 and 64 is to ensure correcttiming of arrival of data within the heuristic processor elements 52,54, 58a and 58b irrespective of data path. The deskewing delay device 66provides for data input synchronously at DI, AI and SAi to give rise tosynchronous outputs at an output interface indicated by a line 68. Theline 68 has portions 68a and 68b corresponding to outputs fromrectangular array portions 58a and 58b respectively. The use oftriangular arrays of clocked latches such as 62, 64 and 66 is well knownin the art of systolic array processing and will not be describedfurther.

The heuristic processor 44 has two alternative modes of operationselectable by a signal applied to an input Tr/Te. A high (binary 1)signal implements Tr (training) mode, and a low (binary 0) signalimplements Te (test) mode. This mode selection signal is passed on totriangular and rectangular arrays 54 and 58 via an appropriate delay(not shown) within the arithmetic unit 52 and via a connection 70. Inthe training mode, training data is derived from the nonlinear system 14via the computer 22 and FIR filter 30 as previously described, and isfed to the data input bus DI after delay at 46. The undelayed version ofthe training data is fed to the answer input AI, the switch 42 being setto connect buses 40a and AI together. The delayed data undergoes atransformation in the digital arithmetic unit 52. This involvescomputation of the Euclidean norm (geometrical distance) of each datavector from each of a set of centres (geometrical origins). Here eachdata vector is a respective set of phase space coordinates. Each suchnorm then undergoes a nonlinear transformation. The results of thetransformation then undergo QR decomposition in the triangular array 54.This provides rotation parameters for application to undelayed data(leading by p clock cycles) in the rectangular array 58 received via thedelay device 60. In the training mode, the boundary and internal cellsB₁₁ to I_(r),r+d+1 all operate adaptively; ie in generating and applyingrotation parameters respectively they update their stored matrixelements on every clock cycle in accordance with each new set of inputs.

The result of operation in the training mode is that data delayed by pclock cycles (interval pT) becomes transformed and least squares fittedto undelayed data. The accuracy of the fit is indicated by residualvalues (errors) appearing at the heuristic processor output 68. Anacceptably accurate fit is indicated by small fluctuating residualvalues.

At the end of training, the boundary and internal cells B₁₁ toI_(r),r+d+1 have respective matrix elements stored therein. In the first(multicolumn) rectangular array portion 58a, the elements map a pointhaving co-ordinates gk(i) (k=1 to d) on a trajectory in a d-dimensionalphase space on to a point p steps later on that trajectory and havingco-ordinates g_(k) (i+p). In the second (single column) rectangulararray portion 58b, the matrix elements map the point g_(k) (i) (k=1 tod) on to a digitised time series value V(t_(i+n-1+p)) also p stepsahead. When training is complete, the mode control input Tr/Te is set toa low voltage signal level (binary 0) to implement test mode. This isrelayed to the triangular and rectangular arrays 54 and 58 via theconnection 70. It switches all the boundary and internal cells of thekinds B and I in these arrays into modes where internal matrix elementupdate is suppressed. Evaluation and application of rotation parametersin the cells B and I respectively continues however. In consequence, themathematical model expressed by stored matrix elements derived duringtraining is frozen.

There are two alternative test modes of operation depending on thesetting of the switch 42. In both modes, the "known" or "comparison"nonlinear system 14 is replaced by an "unknown" nonlinear systemrequired to be compared with the former. The first mode involves theswitch 42 remaining in the training setting, and thereby connecting busspur 40a to answer input AI. Data derived from the unknown system isthen transformed into co-ordinates g_(k) '(i) etc by the FIR filter 34employing singular vectors f₁ to f_(d) derived from the known system 14.

The co-ordinates g_(k') (i) etc are fed directly to the answer input AIbecause of the setting of the switch 42. They are also fed via the delay46 to the data input DI. The associated time series values V'(t_(i)) arefed to the subsidiary answer input SAI. The input at DI becomestransformed nonlinearly into an r-dimensional space by the arithmeticunit 52. The output from the unit 52 is further transformed by thefrozen mathematical model stored in the triangular and rectangulararrays 54 and 58. After this further transformation it is subtracted inthe rectangular array 58 from co-ordinates and time series valuesreceived from the delay device 60. This gives rise to error values atthe processor output 68. Each error value is the difference between arespective measured value from the unknown nonlinear system and anassociated predicted value calculated by the heuristic processor 44 fromits frozen mathematical model. If the error values are sufficientlysmall the unknown non-linear system is demonstrated to be equivalent tothe known system 14, in that its future behaviour is predictable p stepsahead on the basis of the behaviour of the known system. Thispredictability will occur in both phase space (rectangular array portion58a) and in the corresponding time series (rectangular array portion58b). This is a demonstration that signals from a trial system can beused to detect and quantify their similarity to those from a previouslyobserved "standard" system.

This first test mode of operation is also important for the purpose ofdetecting changes in the operation of systems which maybe chaotic butare deterministic in the mathematical sense. For example, an oscillatorysystem such as a rumbling mechanical bearing may affect readings from atransducer arranged to monitor a parameter of the device of which thebearing is part. In this case, the "known" nonlinear system 14 used intraining will be the device with the bearing in good working order. Thetransducer output (which contains bearing rumble) is also the "unknown"system. The "error values" output at 68 during test mode will be smallso long as the bearing operates normally. The error values will increaseas the bearing degenerates through wear. This is a demonstration thatsignals from systems are separable in accordance with the invention whenthey become dissimilar.

The second test mode of operation requires the switch 42 to be set toits second position. This disconnects the bus spur 40a and applies abinary zero signal level to each segment of the answer input AI. Inother respects, operation in the second test mode is as in the first. Inconsequence, data input at DI and transformed at 52, 54 and 58 becomescompared to zero signals input to the rectangular array 58 from thedelay device 60. In consequence, the error values appearing at theoutput interface 68 are in fact generated from comparisons with zerosignals. These values are therefore estimates of the co-ordinates inphase space (output 68a) and the time series (output 68b) of the"unknown" system p steps (clock cycles) in the future in each case. Thepredictions are made on the basis of the multidimensional matrix modelbuilt up during training of the triangular and rectangular arrays 54 and58 and stored therein. This demonstrates that the prediction device 10has the capacity for generating predictions of future states of anonlinear dynamical system.

The mode of operation of the computer 22 will now be described in moredetail. The computer 22 is programmed to carry out a sequence ofoperations; this produces a singular value decomposition of the timeseries data of the kind V(t_(i)) from the "known" system 14. Thesequence is as follows:

(a) Determine an acceptable value for n, the number of dimensions inphase space;

(b) Carry out a singular value decomposition of time series V(t_(i)) etcto define a basis for the phase space;

(c) Determine noise-corrupted singular functions generated in (b) andomit them to produce a d-dimensional sub-space;

(d) Load singular vectors into FIR filter 34;

(e) Activate FIR filter 34 by clock signal to read in V(t_(i)) etc; and

(f) Activate heuristic processor 44 n clock cycles after (e) to startinput of co-ordinates and signals of the kind g_(k) (i) and V(t_(i)).

To carry out the above operations, the computer 22 reads in thedigitised time series having the general term V(t_(i)), where i=1 to M.A typical value for M would be 50,000. It converts the terms in the timeseries to a set of normalised values by calculating their statisticalmean and variance, subtracting the mean from each term and dividing theresult by the variance. For convenience, a normalised time series valuederived from and equivalent to a sampled ADC output signal V(t_(i)) willalso be referred to below as V(t_(i)).

For test purposes, a "known" nonlinear system 14 was constructed. Thisconsisted of a mechanical bearing in the form of a ball race, analuminium bar contacting the centre of the race, and a positionaltransducer contacting the bar. The transducer was a piezoelectriccrystal pickup of the kind used to contact records in domestic recordplayers. It gave an electrical output indicating bearing rumble, ieunwanted and at least. partly deterministic nonlinear bearing vibration.

The transducer produced a time series of 50,000 (ie M) values. The first2048 of these are shown in FIG. 2 in digitised form after normalisationby the computer 22. To carry out operation (a) (as defined above) anddetermine an acceptable number n of phase space dimensions, the computer22 calculates the autocorrelation function of the (now normalised) timeseries V(t_(i)) etc. Writing V_(i) for V(t_(i)), the general term of thenormalised time series, in order to reduce complexity of equations, theautocorrelation function G has successive values of the kind G_(k)defined by: ##EQU2## where k=0, 1, 2, 3, . . . (K-1), where K is themaximum value of k, and K is less than M.

FIG. 3 shows the autocorrelation function G for the normalised timeseries data shown in FIG. 2. The computer 22 calculated the function Gover values of k from 1 to 128, and the values G_(k) are shown plottedin arbitrary units against k in FIG. 3.

The autocorrelation function G crosses the k axis, ie G_(k) =0, atapproximately k=28. It has been found empirically that an appropriatevalue of n, the number of phase space dimensions, is equal to the firstvalue of k at which G_(k) =0. The value of n is therefore chosen to be28.

The computer 22 now implements operation (b), the singular valuedecomposition of the time series. To do this, it forms a matrix X inwhich successive rows are successive Takens' vectors T_(i) etc. Thelatter are defined by

    T.sub.i =[V.sub.i, V.sub.i+1, . . . V.sub.i+n-1 ]          (5)

Equation (5) shows that the ith Takens' vector T_(i) (i=1 to M-n+1) hasvector elements which are a sequence of n successive time series valuesbeginning with V_(i). This is discussed by F Takens in "DetectingStrange Attractors in Turbulence" Lecture Notes in Mathematics Ed. D ARand and L. -S. Young, Springer, Berlin, 1981 p.366. X is thereforegiven by: ##EQU3## The transpose of X is X^(T) given by: ##EQU4## Theautocovariance matrix Ξ is defined by:

    Ξ=X.sup.T X                                             (9)

Ξis then processed by the computer 22 to provide the required singularvalue decomposition. This matrix is subjected to the well knownmathematical process of diagonalisation. The process consists ofHouseholder transformation of the matrix to tri-diagonal form, followedby application of the QR algorithm to provide eigenvalues. Theeigenvalues are the squares of the required singular values. This isdescribed by W H Press et al. in "Numerical Recipes--The Art ofScientific Computing", Cambridge, p.350 et sequi. Also, thedecomposition yields a set of n (in the present example n=28) singularvectors expressed as sets of digital values defining mathematicalfunctions f₁ to f_(n). Each singular vector is associated with arespective singular value s₁ to s_(n), which might loosely be referredto as the strength of the respective vector. The decomposition may bethought of as generating the higher dimensional equivalent of anellipsoid having axes directed as indicated by respective singularvectors and axial lengths given by respective singular values.

This is discussed by W H Press et al (as referred to above) at pages 52to 64.

The lengths or strengths of the singular vectors, ie their respectivesingular values, are arranged to reduce monotonically as j increases,where j is the vector index number in f_(j) and j=1 to n. It willnormally fall to the nonlinear system's noise level before j reaches n.This indicates that singular vectors associated with singular valueswhich are small are significantly corrupted by noise. Such vectorsshould therefore be omitted. Consequently, the dimensionality of thephase space reduced from n to d, where singular values s₁ to s_(d) aregreater in magnitude than the noise level.

In a correctly designed analyser 10, the ADC 16 would be arranged torender the noise level approximately equal to but greater than the leastsignificant bit of the ADC's output digital word. The programming of acomputer to deal with equations (4) to (9), to truncate the singularvalue decomposition and to supply appropriately timed data and clocksignals is straightforward and will not be described. Singular valuedecomposition in the foregoing context is described by Broomhead et alin Physica 20D, 1986, p217.

Referring now to FIG. 4(a) and 4(b), there are shown the results of asingular value decomposition of the 50,000, normalised digital timeseries values obtained from the ball race bearing referred to earlier. Aplot indicated generally by 80 provides the spectrum of the singularvalues s₁ to s_(n) produced in the decomposition. In the present examplen=28, so the values are s₁ to s₂₈ with a general value s_(i) say, inwhich i=1 to 28. The plot 80 shows a monotonically falling curve 82 oflog₁₀ s_(i) against singular value (or singular vector) number i from 1to 28. At i=1, the curve 82 has a logarithm of 0.6 approximately,corresponding to a relative

magnitude of about 4. At i=28, it also falls to a logarithmic value ofabout -2.3, equivalent to about 5×10⁻³. It does not cross the noiselevel of the ADC 16 (indicated by a chain line 84). However, asindicated at a point 85 where horizontal and vertical chain lines 86 and87, an elbow appears in the curve 82. To the left of the vertical line87, the curve 82 varies rapidly with singular value number i; to theright of this line, the variation of log₁₀ s_(i) with i isinsignificant. The vertical line 87 is at i=7. For i equal to 8 or more,it can be inferred that the nineteen singular values are corrupted bynoise arising in the system 14 and sensor transducer 12. This inferencecan be taken because noise is homogeneously distributed, and lack ofvariation in a spectrum is normally associated with a noise level. Thehorizontal chain line 86 is consequently treated as the system noiselevel (as opposed to the ADC quantisation noise level 84). Singularvectors for which i>7 are therefore considered to be unacceptablycorrupted by noise. Singular vectors or functions f₁ to f₇ have anacceptable degree of noise. These are consequently retained, and f₈ tof₂₈ are discarded. This creates a noise-reduced subspace of thenonlinear system, which subspace in the present example has sevendimensions. The value of the parameter d in equation (2.d) is therefore7.

FIG. 49a) and 4(b) also shows the twenty-eight singular vectorsgenerated by singular value decomposition in the present example, thevectors being shown as curves labelled f₁ to f₂₈.

To compare FIG. 49a) and 4(b) with classical Fourier analysis, it isnoted that the latter assumes infinitely long (and thereforeimpractical) time series. Fourier analysis characterises linear systemswith discrete frequencies (line spectra). It is not useful for nonlinearsystems which exhibit broad band spectra.

If required, it is possible to combine the alternatives provided by theswitch 42 in FIG. 1. This requires the switch 42 to be removed and athird portion 58 (not shown) robe added to the rectangular array 58. Theportion 58c would be of like construction to portions 58a and 58bcombined. In training mode, the portion 58c would receive input from AI,which would be bifurcated for this purpose. In test mode the portion 58cwould receive zero inputs. A third output (say 68c) would then providepredicted values of co-ordinates and time series while errors betweenthese values and measured parameters would appear at 68a and 68b.

Referring now to FIG. 5(a) and 5(b ), there are shown results similar tothose shown in FIG. 4 but obtained from a modified Van der Polelectronic oscillator. Such an oscillator has very small inherent noise.Any noise which appears in the course of signal processing is due to ADCquantisation. FIG. 5(a) and 59b) shows a singular value spectrum 90falling below an estimate of an ADC quantisation level 91 at i=13approximately. The level 91 is the standard deviation of the relevanttime series divided by one in the associated ADC's least significant bitchannel. The spectrum 90 levels out at about i=20. ADC quantisationnoise therefore dominates the results, and the spectrum requirestruncation at i=13 (giving d=13). With a more sensitive ADC 16,truncation could be extended to i=20.

FIG. 5(a) and 5(b) also show the forms of the singular vectors orfunctions obtained in singular value decomposition. There arethirty-eight, i.e. i=1 to 38.

Referring now to FIG. 6, there is shown a functional block diagram ofpart of a further embodiment t10 of the invention. It is a dynamicalsystem analyser which is a modification of that shown in FIG. 1, andparts equivalent to those described earlier have like references. In thecase of numerical references, 100 has been added to each.

In view of the similarities between the analysers 10 and 110, thedescription of the latter will concentrate on aspects of difference.

The analyser 110 has an heuristic processor 144 and FIR filter 134equivalent to those previously described. The FIR filter 134 isillustrated laterally inverted such that its latch chain 136 appears onthe right. The latch chain 136 is extended by three additional latches136_(n+1), 136_(n+2) and 136_(n+3) as compared to latch chain 36. Thelatches 136_(n) to 136_(n+3) are connected to respective inputs SAI₀ toSAI₃ collectively forming a subsidiary answer input SAI. The inputsSAI₁, SAI₂ and SAI₃ are connected via respective input delay units 115₁to 115₃ to a delay device 160 incorporated in the heuristic processor144. There is no equivalent of the switch 42. The input delay units115₁, 115₂ and 115₃ implement respective delays of τ, 2τ and 3τ, where τis a clock cycle as before.

Signals pass from the delay device 160 into respective columns C₀ to C₃forming in combination an output array 158b. The array 158b isequivalent to a multi-column version of the column 58b of the earlierembodiment, in which least squares fitting to normalised time seriesdata was implemented. The columns C₀, C₁ and C₂ are in series withoutput delay units 1170, 1171 and 1172 respectively; the latter providedelays of 3τ, 2τ and τ respectively, and in consequence signals reachingan output interface 168 are temporally deskewed (as compared to theequivalent output from columns C₀ to C₃). The system analyser 110operates as follows. A training procedure is employed equivalent to thatdescribed earlier. A comparison system (not shown) provides a timeseries of signals from which singular vectors f₁ to f_(d) are derived bysingular value decomposition and truncation. The vectors are loaded intothe FIR filter 134, and the time series V₁ . . . V_(i) . . . is clockedalong the latch chain 136. In training mode, phase space coordinatesoutput from the FIR filter 134 become nonlinearly transformed, QRdecomposed and fitted to time series data as described for the earlierembodiment 10. However, in the present embodiment, instead of a singlenumber (p) of clock cycle delays in the device 46, there are fourindependent delays of 0, 1, 2 and 3 clock cycles (intervals τ). Toimplement this, the delay circuitry has been rearranged somewhat, makinguse of the latch chain 136 and adjusting for delays experienced byrotation parameters passing between successive columns C₀ to C₃. Theinput delay units 115₁ etc provide that adjustment.

In consequence, during training the ith set of phase space co-ordinatesg₁ (i) to g_(d) (i) output from the FIR filter 134 becomes associated bythe heuristic processor 144 with each of the following time seriespoints: V_(i+n-1), V_(i+n), V_(i+n+1) and V_(i+n+2). This is carried outin respective columns C₀ to C₃ ; ie column c_(q+1) is responsible forfitting the transformed version of g₁ (i) to g_(d) (i) to V_(i+n+q),where q is -1, 0, 1 or 2 and q is equivalent to p-1.

The columns C₀ to C₃ provide error values between transformed versionsof g₁ (i) to g_(d) (i) and V_(i+n-1), V_(i+n), V_(i+n+1) and V_(i+n+2).The prediction interval increases from 0 to 3 clock cycles betweencolumns C₀ to C₃, and consequently the disconformity or error is likelyto increase also.

FIG. 7 is an indication of what might be expected as an output from ananalyser similar to the analyser 110 but having eleven output columns C₀to C₁₀ (not shown). It is a plot 200 of the logarithm to base 10 of theRMS prediction error against prediction interval, the latter being thenumber p of clock cycle intervals (τ); here p goes from 0 to 10, and isthe suffix of the respective output column C₀, C₁, C₂ . . . in eachcase. A graph such as the graph 200 may be employed to identify thedegree of predictability of a nonlinear dynamical system. The RMSprediction error shown in FIG. 7 is seen to increase with increasing pup to some maximum level at which it saturates. The saturation level iscaused by the fact that even if the analyser made the worst guesspossible the prediction error will be limited to some maximum,determined by the training and processing `attractor sizes`.

The result shown in FIG. 7 is calculated as follows. The sequence ofvalues x_(t), t=0 to 4000 are generated by iterating the differenceequation

    x.sub.t =1-α(x.sub.t-1).sup.2 +βx.sub.t-2       (10)

The parameters α and β are constants with values α=1.4, β=0.3. Thesevalues are chosen so that the sequence is chaotic, and a small startingvalue is used for x_(o). The sequence of values of x_(t), where t=1001to 2000, are used to train the dynamical system analyser using 36centers uniformly distributed across a plane bounded by the points (2,2),(2, -2), (-2, -2) and (-2, 2). The sequence of values of x_(t), wheret=3001 to 4000, are then processed, and Log₁₀ of the RMS predictionerror is calculated using the respective set of 1000 error values whichappear from each of the output columns C₀ to C₁₀.

The invention may be implemented in an alternative manner to thatpreviously described. The analyser 10 for example was operated insuccessive training and test modes. It is however possible for atraining procedure on one analyser to produce results employed by adifferent device. Referring to FIGS. 1 and 1a once more, the trainingprocedure in the heuristic processor 44 generates a set of matrixelements. The matrix element values become finalised when training isended, and they are stored on individual internal and boundary cells Iand B. In test mode, these values are frozen, and they implement amathematical (matrix) transformation on data output from the digitalarithmetic unit 52, which itself implements a nonlinear expansion abouta set of centres. The matrix transformation is equivalent tomultiplication by a weight vector or weight matrix; individual elementsof the weight vector or matrix are obtainable from the trained heuristicprocessor. The procedure is set out in PCT App. No. PCT/GB90/00142previously cited, and will not be described in detail. However, havingsuch weighting elements and singular vectors determined in a trainingprocedure, together with the centres and nonlinear transformation (unit52) used in training, operation in test mode may be by means of a simpleweighting device operating on the output of the digital arithmetic unit52. This enables the apparatus requirements of test devices to bereduced.

It has been found that the invention might, in certain circumstances,prove to be too sensitive to minor departures of a dynamical system fromexpected behaviour. In other words a small departure might give a largeerror. In this case it is possible to employ the invention in adifferent manner, as will now be described.

As has been said, the heuristic processor 44 operating in training modebuilds up stored elements which are related to a weight matrix. Elementsof the weight matrix may be obtained at the output 68 by input of unitvectors to the triangular systolic array 54 as described inPCT/GB90/00142. This enables a dynamical system to be characterised bythe elements of the weight matrix to which it gives rise in theheuristic processor 44. Data from a comparison dynamical system maytherefore be input to an analyser of the invention and employed thereinto generate weight matrix elements for subsequent read out. Theseelements are then available for comparison with weight matrix elementsof any trial system. If the weight matrices of the comparison system andthe trial system are similar, then the two systems are alike.

It is not in fact essential to convert elements stored on internal andboundary cells I and B to weight matrices in order to compare the trialand comparison systems. Since each set of stored elements is directlyrelated to a weight matrix, and the relationship is invariant, storedelements from the trial and comparison systems may be compared with oneanother directly without intervening transformation to weight matrices.

There are two approaches to obtaining the weight matrix of the trialsystem. One approach is to retrain the analyser 10 using data from thetrial system; ie after the analyser 10 has been trained using comparisonsystem data, the comparison system weight matrix is read out and thewhole training procedure previously described is carried out once more.The computer 20 consequently carries out a singular value decompositionof data derived from the trial system, and it loads the FIR filter 34with singular vectors derived in that decomposition. The heuristicprocessor 44 subsequently processes the output of the FIR filter 34 aspreviously described to build up elements stored on internal andboundary cells.

The second approach is for the heuristic processor 44 to operate aspreviously described for the training mode, but the FIR filter 34 isloaded with singular functions derived from the singular valuedecomposition of data from the comparison system. There is no suchdecomposition of trial system data. In consequence, the trial systemdata is transformed to phase space coordinates using comparison systemsingular functions. Subsequently, the phase space coordinates soproduced undergo processing in the heuristic processor as previouslydescribed to build up a trial system set of stored element related to aweight matrix.

We claim:
 1. A dynamical system analyser comprising:(1) means forderiving a time series of signals from a dynamical system, (2) means forgenerating from the time series a set of singular vectors each vectorassociated with a respective singular value which is greater than anoise level of the dynamical system, said set defining a phase space andcorresponding to a subset of a singular value decomposition of the timeseries, (3) means for transforming the time series into a sequence ofsets of phase space co-ordinates by predicting, for each set, arespective Takens' vector of signals from the time series onto the setof singular vectors, and (4) heuristic processing means operating in atraining mode for carrying out a transformation in which said sequenceundergoes QR decomposition and least squares fitting to reference datain order to generate a mathematical model related to thattransformation.
 2. An analyser according to claim 1 wherein theheuristic processing means operates in a test mode for transforming afurther sequence of phase space co-ordinates, and for providingpredictive output values on the basis of the mathematical model.
 3. Ananalyser according to claim 1 wherein the heuristic processing means intraining mode employs reference data comprising at least one of the timeseries and the sequence of sets of phase space co-ordinates derivedtherefrom, the analyser including reference data input means forsupplying such data to the heuristic processing means with a relativelead time interval over input of the said sequence for QR decompositionand least squares fitting.
 4. An analyser according to claim 3 whereinthe heuristic processing means is responsive to input reference data intest mode, the reference data consisting of at least one of:(a) phasespace co-ordinates derived from a trial system and having said lead timeinterval, (b) a time series derived from a trial system and having saidlead time interval, and (c) zero data.
 5. An analyser according to claim1 wherein the generating means is a digital computer programmed to carryout a singular value decomposition of the time series of the dynamicalsystem.
 6. An analyser according to claim 1 wherein the transformingmeans is a finite impulse response filter with an input of the timeseries in the form of a succession of Takens' vectors, and to computethe projections of each Takens' vector on to each of the set of singularvectors for providing a respective set of co-ordinates in a phase spacefor each Takens' vector.